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ABSTRACT 

We describe a methodology to accurately compute halo mass functions, progenitor 
mass functions, merger rates and merger trees in non-cold dark matter universes using 
a self-consistent treatment of the generalized extended Press-Schechter formalism. Our 
approach permits rapid exploration of the subhalo population of galactic halos in 
dark matter models with a variety of different particle properties or universes with 
rolling, truncated, or more complicated power spectra. We make detailed comparisons 
of analytically derived mass functions and merger histories with recent warm dark 
matter cosmological N-body simulations, and find excellent agreement. We show that, 
once the accretion of smoothly distributed matter is accounted for, coarse-grained 
statistics such as the mass accretion history of halos can be almost indistinguishable 
between cold and warm dark matter cases. However, the halo mass function and 
progenitor mass functions differ significantly, with the warm dark matter cases being 
strongly suppressed below the free-streaming scale of the dark matter. We demonstrate 
the importance of using the correct solution for the excursion set barrier first-crossing 
distribution in warm dark matter-if the solution for a flat barrier is used instead the 
truncation of the halo mass function is much slower, leading to an overestimate of the 
number of low mass halos. 
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1 INTRODUCTION 

The cold dark matter (CDM) paradigm (Blumenthal et al. 
1984) works extremely well on large scales (Seljak et al. 2005; 
Percival et al. 2007; Ferramacho et al. 2009; Sanchez et al. 
2009; Komatsu et al. 2011), but there remains the possibility 
of deviations from the CDM expectations on small scales- 
arising notably from the issue of cores vs. cusps and the inner 
slope of dark matter density profiles (Salucci 2001; Donato 
et al. 2004; Newman et al. 2009; Donato et al. 2009; de Blok 
2010; Kuzio de Naray & Kaufmann 2011; Kuzio de Naray & 
Spekkens 2011; Newman et al. 2011; Salucci et al. 2012; Wolf 
& Bullock 2012) and the apparent paucity of bright satel- 
lites around the Milky Way (Boylan-Kolchin et al. 2011, 
2012). Several extensions to dark matter phenomenology 
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(Markevitch et al. 2004; Boehm & Schaeffer 2005; Arm & 
Shapiro 2005; Miranda & Maccio 2007; Randall et al. 2008; 
Boyarsky et al. 2009; Lovell et al. 2012) and several differ- 
ent particle physics candidates for dark matter (Raffelt 1990; 
Turner 1990; Jungman et al. 1996; Hogan & Dalcanton 2000; 
Spergel & Steinhardt 2000; Abazajian et al. 2001; Cheng 
et al. 2002; Feng et al. 2003; Sigurdson & Kamionkowski 
2004; Hubisz & Meade 2005; Feng et al. 2009) have been 
put forward to explain these deviations, although it remains 
unclear if any of these proposals is able to fully explain ob- 
served phenomena (Kuzio de Naray et al. 2010) or if they are 
even necessary, with the observed phenomena simply being a 
consequence of galaxy formation physics in a CDM universe 
(Benson et al. 2002a; Libeskind et al. 2007; Li et al. 2009; 
Stringer et al. 2010; Font et al. 2011; Oh et al. 2011; Gover- 
nato et al. 2012; Kuhlen et al. 2012; Pfrommer et al. 2012; 
Pontzen & Governato 2012; Starkenburg et al. 2012; Vera- 
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Ciro et al. 2012). A variety of experimental measurements 
may be sensitive to the small-scale structure of dark matter 
halos (Simon et al. 2005; Viel et al. 2008). The most promis- 
ing are future lensing experiments, which have the potential 
to strongly constrain dark matter particle phenomenology 
(Keeton & Moustakas 2009; Vegetti & Koopmans 2009a,b; 
Vegetti et al. 2010a,b, 2012). A variety of previous work has 
shown that the subhalo mass function should depend sen- 
sitively on the dark matter physics and the shape of the 
power spectrum. To maximize the scientific return of future 
experiments therefore requires the ability to accurately and 
rapidly predict the distribution of dark matter substructure 
as a function of dark matter particle thermal or interaction 
properties, for arbitrary power spectra. 

In this work, we develop techniques to follow the growth 
of nonlinear structure in non-CDM universes using a fully 
consistent treatment of the extended Press-Schechter for- 
malism. We demonstrate the performance of these tech- 
niques by applying them to a representative case of warm 
dark matter (WDM), which has the advantage of several 
pre-existing N-body simulations which we utilize to test the 
accuracy of our methods. WDM particles are lighter than 
their CDM counterparts, allowing them to remain relativis- 
tic for longer in the early universe and to retain a non- 
negligible thermal velocity dispersion. This velocity disper- 
sion allows them to free-stream out of density perturbations 
and so suppresses the growth of structure on small scales 
(Bond & Szalay 1983; Bardeen et al. 1986). While mass func- 
tions have been previously considered in this case (Barkana 
et al. 2001), we go one step further and develop methods to 
compute conditional mass functions, and halo merger rates 
and use these to construct merger trees in WDM universes. 
These merger trees are a key ingredient required to predict 
the distribution of substructure masses, positions and inter- 
nal structure as is necessary to make detailed predictions for 
future lensing experiments. 

The remainder of this paper is organized as follows. In 
§2 we describe the changes that we introduce to the ex- 
tended Press-Schechter formalism to make it applicable to 
the case of non-CDM scenarios (including some specifics for 
the WDM case). In §3 we apply these methods to the case 
of WDM, first comparing their predictions to the available 
data from N-body simulations, then exploring the limita- 
tions of ignoring the effects of WDM velocity dispersion, 
and presenting a comparison of key results between WDM 
and CDM. Finally, in §4 we discuss the consequences of this 
work and present our conclusions. 

We also include two appendices. Appendix A gives a de- 
tailed derivation of our numerical procedure for solving the 
excursion set first crossing problem for arbitrary barriers. 
Appendix B explores the numerical accuracy and robust- 
ness of the methods developed in this work. 

When comparing our analytic theory with results from 
N-body simulations we will adopt the same cosmological pa- 
rameters and dark matter particle properties as were used 



for the simulation. These values will be listed where relevant. 
For the rest of this work, specifically in §2, §3.2, §3.3 and 
Appendix B we adopt a canonical a cosmological model with 
Q M = 0.2725, fi A = 0.7275, n b = 0.0455 and H = 70.2 km 
s _1 Mpc~ x (Komatsu et al. 2011) and a canonical WDM 
particle of mass, m x = 1.5 keV and effective number of de- 
grees of freedom <?x = 1.5 (the expected value for a fermionic 
spin- 1 particle). 



2 METHODS 

Our approach makes use of the Press-Schechter formalism 
(Press & Schechter 1974; Bond et al. 1991; Bower 1991; 
Lacey & Cole 1993) which, after substantial development 
and tuning against N-body simulations, has proven to be 
extremely valuable in understanding the statistical prop- 
erties of dark matter halo growth in CDM universes. The 
Press-Schechter formalism in its modern form is expressed 
in terms of excursion sets-the set of all possible random 
walks in density at a point as the density field is smoothed 
on ever smaller scales. Halo formation corresponds to a ran- 
dom walk making its first crossing of a barrier. The height 
of that barrier is determined from models of the non-linear 
collapse of simple overdensities. 

The Press-Schechter algorithm requires three ingredi- 
ents: 1) the power spectrum of fluctuations in the density 
field (characterized by a(M), the fractional root- variance 
in the linear-theory density field at z = 0); 2) the critical 
threshold in linear-theory corresponding to the gravitational 
collapse of a density perturbation, <5 C ; and 3) a solution for 
the statistics of excursion sets to cross this threshold. We 
will address each of these three ingredients below. 



2.1 Power Spectrum 

We assume a power-law primordial power spectrum with 
n s = 0.961 (Komatsu et al. 2011), and adopt the transfer 
function of Eisenstein & Hu (1999). We include a modifi- 
cation for warm dark matter using the fitting function of 
Bode et al. (2001; as re-expressed by Barkana et al. 2001) 
to impose a cut-off below a specified length scale, A s : 

T(k) -> T(k) [1 + (ek\ s ) 2 "] ~ V '\ (1) 

where e = 0.361, r\ = 5 and v = 1.2 are parameters of the fit- 
ting function. For our canonical WDM particle, the smooth- 
ing scale 2 is A s = 0.124 Mpc (Barkana et al. 2001; eqn. 4) 
corresponding to a mass of M B = 47rpAj /3 = 2.97 x 10 8 M Q . 
The power spectrum is normalized to give the required 
as = 0.807 (Komatsu et al. 2011) when integrated un- 
der a real-space top-hat filter of radius 8/z _1 Mpc (where 
h = i/o/100 km s^ 1 Mpc" 1 ). 



1 The two usual candidates-both lying beyond the standard 
model of particle physics-arc sterile neutrinos (Dodclson & 
Widrow 1994; Shaposhnikov & Tkachcv 2006) and gravitinos (El- 
lis et al. 1984; Moroi ct al. 1993; Kawasaki et al. 1997; Gorbunov 
et al. 2008). 



2 This scale is usually approximated as being equal to the speed 
of the particles at the epoch of matter-radiation equality multi- 
plied by the comoving horizon scale at that time; see Bode et al. 
(2001) for further discussion. 
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2.1.1 Window Function 

To derive the variance, S(M), or, equivalently, the root- 
variance, a(M) = S(M), from the power spectrum a win- 
dow function, W(k\M), must be adopted. Specifically, 



S(M) =a 2 (M) 



i r 



Ank 2 P(k)W 2 (k\M)dk. 



(2) 



In the Bond et al. (1991) re-derivation of the Press 
& Schechter (1974) mass function and extended Press- 
Schechter conditional mass functions one assumes a sharp 
fc-space filtering, 

1 iffc^fcs(M) 
iffe>fe s (M), 



W(k\M) 



(3) 



of the linear density field so that the trajectories of density 
fluctuation, 8, versus mass scale are true Brownian random 
walks. Despite the use of this sharp fc-space filtering in de- 
riving the extended Press-Schechter solutions, it has been 
conventional to adopt the o~(M) that is given by a real-space 
top-hat window function, 

3[sin(fcfl) - kRcos(kR)] 



W(k\M) = 



(4) 



(fci?) 3 

when utilizing those solutions. The reason for this is that 
only in this case is there no ambiguity between the filtering 
scale, R, and the corresponding mass, 

/3M\ 1/3 



R - 



\4-7rp J 



(5) 



(see Lacey & Cole 1993). For CDM power spectra, which are 
essentially pseudo power-laws with a slowly varying slope, 
the shape of the a(M) function varies little with the choice 
of filter function. 

However, in the case of power spectra with a small scale 
cut-off, such as those of WDM models and many other non- 
CDM candidates, the situation is very different. With a top- 
hat window function, a(M) continues to increase with de- 
creasing filter mass even for masses well below the cutoff 
scale. This happens because although there are no new in- 
trinsic small scale modes entering the filter, the longer wave- 
length modes are getting re-weighted as the mass scale of 
the filter increases. In contrast, with a sharp-fc filter, a(M) 
increases monotonically up to the cutoff scale and then be- 
comes constant with the transition being determined by the 
abruptness of the cutoff in the input WDM power spectrum. 
Since the Press-Schechter mass functions depend directly 
on a(M), including a direct dependence on the logarith- 
mic slope of o~(M), its predictions for the case of WDM are 
sensitive to this choice. Hence for non-CDM models it is im- 
portant to revisit the issue of filter choice and one would 
expect that the conventional top-hat choice will lead to an 
overestimate of the number of low mass haloes (an effect 
apparent in the works of Barkana et al. (2001) and Menci 
et al. (2012)). 

We begin by adopting the sharp fc-space window func- 
tion of eqn. (3) when working with truncated power spectra, 
such that the contribution from long-wavelength modes re- 
mains constant. In the case of a top-hat real-space window 
function there is a natural relation between the top-hat ra- 
dius and the smoothing mass scale, as given by eqn. (5). In 
the case of a sharp fc-space window function, the choice of 
k B (M) is less clear. Lacey & Cole (1993) suggest choosing 



W{r = 0) = 1 (where W(r) is the Fourier transform of the 
window function) and then choosing k B (M) such that the 

integral of the mean density, p, under W(r) equals the re- 
quired smoothing mass, M. However, as noted by Lacey & 
Cole (1993), this choice for fc s (M) lacks any strong physical 
motivation. 

Therefore, we advocate a different approach, namely we 
choose 



fc s = a/R, 



(6) 



where R is computed from eqn. (5) and the parameter a is 
chosen such that the turnover in the halo mass function oc- 
curs in the same location as that seen in N-body simulations 
of WDM, as will be shown in §3.1. We find that a = 2.5 is 
required to meet this condition 3 .. This will mean that a(M) 
differs from the usual result on large scales, resulting in a dif- 
ference in the mass function for high mass halos that would 
not be expected. To remedy this, we note that the critical 
overdensity for collapse is usually motivated on the basis of 
spherical top-hat collapse models. Since we are no longer us- 
ing a top-hat filter we choose to allow freedom in the choice 
of the critical overdensity such that the halo mass function 
is unchanged for high masses. This point will be discussed 
in greater detail in §2.2.2. 

We note that fixing the fc s (M) relation in this way intro- 
duces a free parameter, a, into our method and introduces 
a dependence on the N-body simulation to which we will 
compare our results in §3.1.1. However, given that the CDM 
linear power spectrum is close to a power-law on the scales 
that will be of interest for WDM models, and that the mod- 
ification to the transfer function due to WDM scales simply 
with the mass of the dark matter particle we expect that 
the same choice of a should be valid for all WDM particle 
masses of interest. 

Fig. 1 shows the resulting root variance in the density 
field (smoothed using a top-hat filter in real-space) as a func- 
tion of mass enclosed within the filter for both WDM and 
CDM models. As expected, a(M) in our WDM model is 
suppressed below the CDM value for M < M s . 



2.2 Critical Overdensity for Collapse (Excursion 
Set Barrier) 

The Press-Schechter algorithm as formulated by Bond et al. 
(1991) associates the collapse of a dark matter halo with a 
random walk in S(M) making its first crossing of some bar- 
rier, B(S). This barrier is usually associated with the criti- 
cal linear theory overdensity for collapse of spherical top-hat 
perturbations (although see §2.2.1 for a refinement of this 
assumption). In the CDM case, that critical overdensity is 
independent of mass scale (since there are no preferred scales 
in the problem). In the non-CDM case, this is no longer true. 
For example, with WDM we expect collapse to become sig- 
nificantly more difficult on small scales, where the WDM 
particles can stream out of collapsing overdensities due to 



3 The normalization advocated by Lacey & Cole (1993) is equiv- 
alent to a = (97I-/2) 1 ' 3 SB 2.42 which is not too different from our 
best-fit value of 2.5 
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M [M G ] M/M, 



Figure 1. The fractional root variance of the mass density field 
as a function of the mean mass contained within that filter. Both 
CDM and our canonical WDM cases are shown. A top-hat real- 
space window function is used for the CDM while a sharp fc-space 
filter is used for WDM. The WDM case is suppressed below the 
CDM value for M < M s . 



their non-zero random velocities. Barkana et al. (2001) ad- 
dressed the question of collapse thresholds in WDM uni- 
verses by performing a set of 1-D hydrodynamical simula- 
tions, in which pressure acted as a proxy for the velocity 
dispersion of WDM particles. They find that the growth of 
collapsing overdensities is suppressed below a characteristic 
mass scale-i.e. the threshold for collapse increases rapidly 
with decreasing mass below the characteristics mass scale. 

We find that the results of Barkana et al. (2001) can be 
well fit by the functional form 4 : 



5 c ,wdm (M, t) = 5 Cj cDM(i) yi(x) 
+ [1 — h(x)] exp 



0.04 



exp(2.3x) 
0.31687 
exp(0.809a;) 



(7) 



where x = log(M/Mj), M is the mass in question, Mj is the 
effective Jeans mass of the warm dark matter as defined by 
Barkana et al. (2001; their eqn. 10): 

1/2 



Mj 



3.06 x 10° 



1.5 



1 + 2. 



cq 



3000 

m x 



0.15 

M , 



,1.0 keV/ 

the redshift of matter-radiation equality is given by 



Zeq = 3600 



57m 
0.15 



and 



h(x) = 1/{1 + exp[(x + 2.4) /0.1]}. 



(8) 



(9) 



(10) 



4 This fit is accurate for the regime where <5 c ,wdm/<5 c ,cdm < 600, 
but substantially over-predicts the results of Barkana et al. (2001) 
for smaller masses. This is a deliberate choice-our aim was to 
match the shape of the function through the region where it tran- 
sitions away from the CDM value. On smaller mass scales the sup- 
pression is so dramatic that the precise value of <5 c ,wdm/<5c CDM 
is unimportant. 



Figure 2. The critical linear theory overdensity for the collapse 
of a spherical top-hat density perturbation in a WDM universe 
normalized to that in a CDM universe as a function of halo mass, 
M (expressed in units of the effective Jeans mass of the WDM, 
Mj). For M/Mj > 1 the WDM and CDM cases coincide, but for 
M/Mj < 1 WDM requires a much larger overdensity to undergo 
gravitational collapse. 



The ratio of the resulting critical overdensity to the 
CDM value with masses scaled to Mj, is shown in Fig. 2. 
For M/Mj < 1 the critical overdensity for collapse in WDM 
universes is much higher than in CDM as a result of the 
non-zero velocity dispersion of warm dark matter particles. 
A small-scale perturbation must have a much larger density 
for its self-gravity to overcome this velocity dispersion and 
cause collapse. 

We emphasize that the calculations of Barkana et al. 
(2001) are approximate as they are based on a hydrodynam- 
ical approximation which will not fully capture the collision- 
less dynamics of WDM. We expect that the approximation 
made by Barkana et al. (2001) could lead to a small (order 
unity) difference in the characteristic mass scale for sup- 
pression of overdensity collapse (e.g. in a similar way that 
the Toomre stability threshold for collisionless stars and gas 
differ by a factor of 7r/3.36). There could plausibly be dif- 
ferences in detail in the shape of the collapse threshold as 
a function of mass, but we are unable to speculate what 
form those might take. A more realistic calculation using 
a Boltzmann solver should be carried out to improve upon 
these results, and explore the dependence on cosmological 
parameters. 



2.2.1 Barrier Remapping 

In the CDM case the original Press-Schechter algorithm as 
formulated by Bond et al. (1991) adopted a barrier, B(S), 
for excursion sets equal to the critical linear-theory over- 
density for the gravitational collapse spherical top-hat per- 
turbations and which was independent of mass and equal 
to S c (M[S],t) — 8 c ,o/D(t) where D(t) is the linear theory 
growth function 5 and 5 c ,o is the collapse threshold (equal 



5 Note that, as is conventional, we place this time dependence 
into the critical overdensity, such that we can always work with 
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to (3/20)(12tt) 2/3 « 1.686 in an Einstein-de Sitter universe; 
solutions for other cosmologies are given by Percival (2005) 
for example) . It is now well known that this constant barrier 
does not result in a halo mass function that agrees well with 
that measured from N-body simulations of CDM (e.g. Tin- 
ker et al. 2008). Motivated by this discrepancy, Sheth et al. 
(2001) introduced a remapping of the barrier 6 calibrated to 
improve the match: 



B(S) -> B{S)VA ( 1 



g 

AB 2 (S) 



(11) 



where A = 0.707, b = 0.5 and c = 0.6 are parameters that 
were tuned to obtain the best match. We retain this remap- 
ping of the barrier in this work 7 when computing halo mass 
functions but not when computing halo merger rates (see 
§2.4 for further discussion of this point). 



2.2.2 Barrier Scaling 

As noted in §2.1.1 our choice of normalization for the window 
function used to derived a(M) from P(k) results in a(M) for 
the WDM case lying above the CDM case for large masses. 
This will change the halo mass function for high mass halos- 
something which is not expected (i.e. WDM should behave 
just like CDM on sufficiently large scales). Therefore, we in- 
troduce a global re-scaling of the barrier (applied after the 
remapping described in §2.2.1-an important point since that 
remapping is nonlinear), multiplying it by a factor of 1.197. 
This value is chosen to counteract the higher a(M) on large 
scales 8 and ensure that the mass function of high mass halos 
remains unchanged (see §2.1.1 for discussion of this point). 



cr(M) at z = 0. The growth function used here is the usual growth 
function computed for CDM, consistent with the definitions of 
Barkana et al. (2001). 

6 This remapping was motivated by considerations of ellipsoidal 
collapse, the barrier for which differs from that for spherical col- 
lapse. In particular, Sheth et al. (2001) noted that the density 
perturbations leading to low-mass halos are expected to be more 
ellipsoidal than those that give rise to the most massive halos- 
effectively making the collapse barrier a function of mass scale. 

7 This mapping was derived strictly for the CDM case. We retain 
it here since for halos with masses, M, much greater than M B we 
expect the WDM case to converge to the CDM solution. How- 
ever, there is no guarantee that this mapping remains accurate 
for M < M B . In particular, it is possible that low mass halos close 
to the cut-off scale could be more spherical in WDM than in CDM 
due to the isotropizing effects of the WDM velocity dispersion. 
Calibration of our methods against reliable N-body simulations 
of WDM universes would be required to evaluate the accuracy 
of the remapping in this regime. The N-body simulations cur- 
rently available and examined later in this work do not address 
this particular issue as they do not include the effects of the non- 
zero velocity dispersion of WDM. However, it is understood how 
to include velocity dispersions in WDM simulations (Brandbygc 
et al. 2008), and this should be included in future simulations. 

8 That is, the ratio of cr(M) computed using our sharp fc-space 
filter to that computed using the top-hat filter is \/1.197. Note 
that this scaling factor does not constitute a free-parameter in 
our approach. Instead, it is fixed to offset the change in the S(M) 
relation on large scales resulting from the difference between top- 
hat and sharp fc-space window functions. Therefore, given a power 
spectrum and a value of a (the coefficient in eqn. 6), the scaling 
factor is directly computable and uniquely determined. 
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M[M Q ] 

Figure 3. The barrier for excursion sets in both the CDM and 
our canonical WDM cases. For large masses the two are offset by 
the constant factor of 1.197 which accounts for the difference in 
definition (top-hat real space versus sharp /c-space window func- 
tions) of <r 2 (M) in the two cases, but the WDM barrier becomes 
very large for M M s . Below M ~ M s the WDM barrier rises 
rapidly due to the velocity dispersion of WDM particles. 



This re-scaling of the barrier is always included, both when 
computing halo mass functions and also when computing 
halo merger rates. We note that this factor is very insensi- 
tive to cosmological parameters and as as it depends only 
upon the shape of the power spectrum (not the normaliza- 
tion). For example, if we compute the appropriate rescaling 
factor for WMAP-1 cosmological parameters (rather than 
the WMAP-7 cosmological parameters used throughout this 
work) we find that the factor changes by < 0.1% (despite 
there being a difference of 11% in erg). 

The resulting excursion set barriers for both CDM and 
our canonical WDM model (including both remapping and 
subsequent re-scaling) are shown in Fig. 3. For large masses 
the two are offset by the factor of 1.197, but the WDM 
barrier becomes very large for M <C M s due to the velocity 
dispersion of the WDM particles. 



2.3 Excursion Sets/Barrier Crossing 

Given a barrier, B(S), and the variance of the density field, 
S(M), the Press-Schechter algorithm proceeds by following 
random walks in S(S) beginning from (S, S) = (0, 0). When a 
given random walk first crosses the barrier at some variance 
S, it is assumed that the point in question has collapsed into 
a halo of mass M(S). The fraction of random walks crossing 
between S and S + dS, f(S)dS, corresponds to the fraction 
of mass in the universe in halos of mass M to M + AM. 

We therefore must compute the probability for random 
walks to cross our barrier between S and S + dS. We will as- 
sume a Gaussian distribution of density perturbations, mo- 
tivated by the simplest inflation models. In the case of a 
constant (or linear in S) barrier, there is a well-known ana- 
lytic solution (Bond et al. 1991; Lacey & Cole 1993; Sheth 
1998). However, for an arbitrary barrier, the solution must 
be computed numerically. Our approach, described in detail 
in Appendix A, is similar to that of Zhang & Hui (2006) 
but is numerically more robust. Briefly, the excursion set 
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problem requires finding a solution to an integral equation: 

rS rB(S) 

1= / f(S')dS' + / P(5,S)dS, (12) 

JO J -ex 

where P(5, S)d5 is the probability for a trajectory to lie 
between 8 and 8 + d8 at variance S. This equation expresses 
mass conservation, i.e. at any S, all trajectories must either 
have crossed the barrier at some smaller S (the first term 
in eqn. 12), or be below the barrier having never crossed at 
smaller S (the second term in eqn. 12). To numerically solve 
this equation we discretize the variance, S, using a grid that 
is uniform in S with a spacing of AS. We then numerically 
solve eqn. (12) as described in detailed in Appendix A. The 
solution is extended to the maximum value of S, if such 
exists (as is the case in a WDM model) or to sufficiently 
large S that smaller mass scales are of no interest for the 
problem in question. The accuracy of our numerical solver 
is explored in Appendix Bl. 

2.4 Merger Rates and Progenitor Mass Functions 

The original Press-Schechter algorithm has been extended 
to compute conditional mass functions (i.e. the mass func- 
tion of halos which will all belong to a single halo of larger 
mass at some later time; Lacey & Cole 1993). Further, the 
form of the conditional mass function in the limit of small 
timestep has been used to estimate merger rates of dark 
matter halos and so to construct merger trees (Cole et al. 
2000). We therefore wish to compute the conditional mass 
function in the WDM case, specifically in the limit of small 
St. In the excursion set approach, the conditional mass func- 
tion is found by simply solving the first crossing problem 
beginning from the (8, S) of the parent halo and using a 
barrier corresponding to an earlier time. This is equivalent 
to solving the original first crossing problem with a modified 
barrier B'(S",ti,t ) = B(S' + S,h) - B(S,t ). 

To compute merger rates we therefore solve the first 
crossing problem using our numerical method but with an 
effective barrier B'(S', ti, to). We choose t\ — (1 — e)to where 
e«1. The rate of crossing this effective barrier can then be 
estimated as f(S)/eto- We will explore the sensitivity of our 
results to the value of the numerical parameter, e. 

In this case the choice of grid in S" is particularly impor- 
tant. Accurate numerical solution requires many grid points 
at small S" (since f(S') will peak very close to S' = for 
small values of e), but also many points close to the max- 
imum possible value of S'. S(M) becomes almost indepen- 
dent of M close to the maximum value of S-as a result many 
points are required to resolve the cut-off of the halo mass 
function as a function of mass. Therefore, we adopt the fol- 
lowing distribution for Sf. 



So = 

Si>o 



1/Si,nn + 1/rSi, 



(13) 



where S^im and 5i,i°g are se ^ s °f points spaced uniformly in 
S and the logarithm of S between the required minimum and 
maximum values and r is a numerical parameter. For r > 1 
the resulting distribution of points is spaced uniformly in the 
logarithm of S for small S, and transitions to being uniform 
in S for large S-as a result, the grid has good resolution both 
close to S — and close to the maximum value of S. The 



value of r controls the location of the transition between 
these two regimes. We find that a value of r = 10 works 
well. The accuracy of our numerical solver for merger rates 
is explored in Appendix B2. 

When solving for halo merger rates we do not in- 
clude any remapping of the barrier function (as discussed in 
§2.2.1). The remapping function was chosen to result in good 
agreement between the Press-Schechter halo mass function 
and that measured in the Millennium N-body simulations 
(Springel et al. 2005). However, our merger rates will be 
used in merger tree construction algorithms which, in the 
CDM case, have been developed to work with a barrier with 
no remapping, but with their own empirical modification of 
merger rates designed to match progenitor mass functions 
found in N-body simulations (see §2.5 for further discussion). 
Scaling of the barrier (see §2.2.2) is included however. 

2.5 Merger Tree Construction 

To build merger trees, we follow the algorithm of Parkinson 
et al. (2008). Briefly, at each point in a merger tree, this 
algorithm evaluates the probability per unit time of a binary 
merger occurring along the branch 



du> 



J 

J A 



M/2 M df dS | dt 
M' dt dM> du> 



G[oj,a(M),a(M')]dM', 



(14) 

where uj = 8 c x)/D{t) and <5 Cj0 is the critical overdensity for 
collapse in the CDM case and the integration is from the 
lowest mass halo to be resolved in the tree, M m i n , to halos 
of half the mass of the current halo, corresponding to an 
equal mass merger. Cole et al. (2000) discuss the subtleties 
of why the integration is carried out over the lower half of 
the progenitor mass range. The rate of accretion of mass in 
halos below the resolution limit of the merger tree 



cLR 
do; 



/" 

Jo 



df dS 
dt dM' 



dt 
duj 



G[oj,a(M),a(M')]dM', (15) 



where in this case the integral is taken over all unresolved 
halos (i.e. those less massive than M m i n ). In the above, 
G[lo, a(M), a(M')] is an empirical modification of the pro- 
genitor mass function introduced to obtain results consis- 
tent with those from N-body simulations of CDM universes. 
Parkinson et al. (2008) show that the form: 



G[UJ, (Ti, CT 2 ] = Go 



(71 



(16) 



with Go = 0.57, 71 = 0.38 and 72 = —0.01 provides a good 
match to progenitor mass functions measured from the Mil- 
lennium Simulation (Springel et al. 2005), which assumes 
similar cosmological parameter values as are used in this 
work. The validity of this empirical modification on much 
smaller mass scales than are probed by the Millennium Sim- 
ulation is explored in Appendix B3. 

We retain this same empirical modification in this work 
since it should remain valid for WDM in the limit of high- 
mass halos. For masses comparable to and below M B this 
empirical modification may no longer be accurate. Our justi- 
fication for retaining this empirical modification is the good 
agreement achieved with progenitor mass functions from N- 
body simulations of WDM as will be shown in §3.1. 

A timestep is then chosen that is sufficiently small such 
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that the probability of multiple branchings is small and the 
fractional change in mass due to subresolution accretion is 
small. This timestep is then taken, branching if a random 
deviate lies below the branching probability, and with mass 
removed at the rate of subresolution accretion. If branching 
does occur the mass of one of the branched halos is selected 
at random from the distribution d 2 f /duidM' and the other 
is chosen to ensure mass conservation. 



2.5.1 Smooth Accretion 

In the CDM case, all mass is locked into halos on some 
scale 9 . This implies that the progenitor mass function in 
CDM contains the entire mass of the parent halo, i.e.: 



/(S')dS' = 1. 



For the WDM case this is not true and we find, 

»s max 



/(S')dS' < 1, 



(17) 



(18) 



i.e. not all trajectories cross the barrier. We identify these 
trajectories as corresponding to smooth accretion, i.e. ac- 
cretion of mass which does not belong to any collapsed 
halo. This is physically distinct from the unresolved halos 
accounted for by eqn. (15), but we can nevertheless account 
for this smooth accretion by boosting d_R/dcj by an amount 
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§3.1.1: 0.25keV - Schneider et al. (2012) N-body 
§3.3 : 1.50keV - CDM vs. WDM comparison 
§3.1.2: 2.20keV - Aquarius WDM counterparts 
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Figure 4. Dark matter halo mass functions, computed using the 
methods described in §2, for the three different WDM particle 
masses considered in this work. The labels give the particle mass 
and the section of this work in which that mass is used. 



comparisons of WDM and CDM. For reference, Fig. 4 shows 
the halo mass function in these three cases, illustrating the 
expected increase in the cut-off mass, m B , as the WDM par- 
ticle mass is decreased. 



dR 
dw 



dfn 

At 



dt 
dw 



(19) 



where f n = 1 — J Q Smax /(S')dS". Note that we choose to 
include the empirical correction of Parkinson et al. (2008) 
here, using the maximum value of a(M) for the ac argu- 
ment. This ensures a treatment consistent with resolved ha- 
los but, once again, should be tested and calibrated against 
WDM simulations. The importance of this smooth accretion 
is explored in Appendix B4. 



3 RESULTS FOR WARM DARK MATTER 

The methods described in §2 have been implemented 
within the open source semi-analytic galaxy formation 
code, GALACTICUS (Benson 2012). All results presented 
in this section were generated using GALACTICUS vO . 9 . 1 
r903. Control files and scripts to generate all results pre- 
sented in this paper using GALACTICUS can be found at 
http: //users . obs . carnegiescience . edu/abenson/ 
galacticus/parameters/dmMergingBeyondCDMl .tar .bz2. 

We consider three different WDM particle masses. In 
§3.1.1 we use mx = 0.25 keV to match the N-body simu- 
lations of Schneider et al. (2012), in §3.1.2 we use mx = 
2.2 keV to match the Aquarius WDM counterpart simula- 
tions (described below), and in §3.3 we use mx = 1.5 keV in 



9 Assuming that cr(M) continues to rise monotonically to arbi- 
trarily small scales. In practice this is not true as even CDM 
will have some cut-off in its power spectrum on very small scales 
(Green et al. 2004; Loeb & Zaldarriaga 2005). However, for our 
present purposes the assumption that all mass is locked into halos 
in CDM is a good one. 



3.1 Comparison with N-body Simulations 

3.1.1 Halo Mass Function 

N-body simulations should, in principle, provide an accu- 
rate determination of the dark matter halo mass function 
in WDM cosmologies, provided that initial conditions are 
constructed carefully. Points in Fig. 5 show the mass func- 
tion measured in an N-body simulation of 0.25 keV WDM 
carried out by Schneider et al. (2012). The upturn below 
2 x 10 11 Mq is a numerical artifact, arising from the frag- 
mentation of filaments due to particle discreteness (Wang 
& White 2007). This is a challenging problem for N-body 
simulations of WDM as the mass scale at which the upturn 
appears decreases only as A r_1//3 (where N is the particle 
number) . 

The solid lines in Fig. 5 shows the mass function pre- 
dicted by our calculation using the same cosmological and 
WDM particle parameters as Schneider et al. (2012). Schnei- 
der et al. (2012) identified halos using a friends-of-friends 
algorithm with a linking length of b = 0.2, corresponding 
approximately to density contrasts of 200. In our model, the 
relevant density contrasts are those arising from the spher- 
ical top hat collapse model (e.g. Percival 2005). We correct 
the masses reported by Schneider et al. (2012) for this dif- 
ference by assuming the halos have NFW density profiles 
(Navarro et al. 1997) with concentrations given by the CDM 
results of Gao et al. (2008) modified by the WDM-to-CDM 
conversion factor reported by Schneider et al. (2012). Addi- 
tionally, Schneider et al. (2012) do not include any velocity 
dispersion of dark matter particles in their initial conditions 
and so the fairest comparison is one in which we do not 
modify the critical overdensity for collapse as described in 
§2.2. 

Results of three different calculations from this work 
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Figure 5. The dark matter halo mass function for WDM in the 
case of a 0.25 keV particle. Points show the results of an N-body 
simulation of WDM from Schneider et al. (2012) (small points are 
used to indicate the region in which the N-body simulation re- 
sults are unreliable as a result of being dominated by halos formed 
through artificial fragmentation), while the line shows the result 
from this work with dark matter properties and cosmological pa- 
rameters matched to those used by Schneider et al. (2012). Note 
that the upturn in the N-body mass function below 2 X 10 Mq 
arises to artificial fragmentation of filaments (see Wang & White 
2007). 

are shown in Fig. 5. Going from top to bottom, the first two 
lines do not include the effects of velocity dispersion (consis- 
tent with Schneider et al. 2012). The first line corresponds 
to a case where we use a top-hat real-space window function 
to determine a(M). It can be seen that this curve, while in 
good agreement with the N-body results for high masses, 
does not produce sufficient suppression of the mass function 
at lower masses-a similar discrepancy was seen by Barkana 
et al. (2001) when comparing their Press-Schechter-based 
model for WDM halo formation to the N-body simulations 
of Bode et al. (2001). The second line (blue) switches to us- 
ing a sharp fc-space window function as described in §2.1.1. 
For high mass halos, this results in only a small change in the 
mass function compared to the CDM case due to the differ- 
ence in S(M) when computed with sharp fc-space and top- 
hat window functions (even after scaling the barrier height 
to compensate for this difference as much as possible; see 
§2.2.2) 10 , but results in an excellent match to the suppres- 
sion of the abundance of low mass halos, indicating that the 
discrepancies found in previous works were due to the artifi- 

10 The slight difference between this curve and that computed 
using a top-hat window function is due our inclusion of the Sheth 
et al. (2001) remapping. If this were not included, the factor of 
1.197 increase the barrier would precisely compensate for the dif- 
ference in o(M) between these two curves on large scales. With 
the Sheth et al. (2001) remapping included this constant factor 
cannot correct for the offset fully (due to the non-linear nature of 
the Sheth et al. (2001) remapping). Future, high-precision work 
should consider re-calibrating the Sheth et al. (2001) remapping 
to match a CDM halo mass function with cr(M ) computed using 
our sharp fc-space window function. This would obviate the need 
for a separate multiplicative increase in the barrier and could 
better capture the scale-dependence of the required correction. 



cial increasing in a(M) below the cut-off which arises when 
a top-hat real-space filter is used. At the lowest masses the 
suppression of the mass function is masked in the N-body 
simulation as it is masked by the upturn due to artificial 
fragmentation of filaments. Finally, the green line adds in 
the effects of velocity dispersion (which are not included in 
the N-body simulation), illustrating the importance of this 
effect to accurately model the suppression of the lowest mass 
halos. 



3.1.2 Progenitor Mass Functions/Merger Rates 

A set of WDM Milky Way mass haloes have been simulated 
and analyzed in Lovell et al (in prep) . The haloes they sim- 
ulated are WDM counterparts of the Aquarius CDM haloes 
presented in Springel et al. (2008), and so have masses of 
order 1O 12 M -we will therefore refer to them as "Milky 
Way-mass WDM simulations". Here we make use of a set 
of four (A-D) haloes simulated with approximately 40 mil- 
lion particles within their virial radii (resolution level 3 in 
the notation of the Aquarius project). One of the haloes 
(A) has also been run at higher (level 2) resolution and we 
have used this to verify that the progenitor mass functions 
that we show in Fig. 6 have accurately converged over the 
range of masses plotted. The cosmological parameters for 
this set of simulations were chosen to match the WMAP7 
results of Komatsu et al. (2011), however the CDM power 
spectrum was modified using the prescription of Bode et al. 
(2001) for WDM as given in our eq. (1) with a smoothing 
scale of A s = 78.4 kpc, corresponding to an approximately 
2.2 keV thermal WDM particle. In each of the simulations 
haloes and subhaloes were identified at each of the 128 out- 
put times using the FOF (Davis et al. 1985) and SUBFIND 
(Springel et al. 2001) algorithms respectively. Merger trees 
of subhaloes were constructed by identifying the descendant 
of each subhalo and halo merger trees built from the sub- 
halo membership of each halo as described in Jiang et al. 
(in prep) and Merson et al. (2012). These halo merger trees 
can be thought of as FOF halo merger trees that have been 
cleaned up to avoid problems that occur when FOF haloes 
are essentially composed of two distinct haloes linked by a 
low density bridge. 

The issue of spurious haloes formed by the fragmen- 
tation of filaments due to numerical noise (Wang & White 
2007) is investigated in Lovell et al. (in prep). They find 
that most of such (sub)haloes can be identified by looking 
at the shape in the initial conditions of the region from which 
their particles originated. The spurious (sub)haloes typically 
originate from very flattened configurations. Lovell et al. de- 
termine a threshold on the axis ratio, c/a, of the inertia ten- 
sor of the initial particle distribution such that for a CDM 
simulation 99% of haloes pass this cut while in a WDM sim- 
ulation most of the spurious haloes are excluded. Here we 
use this criterion to exclude complete subhalo branches from 
the N-body merger trees. In these N-body merger trees the 
same subhalo is identified at subsequent epochs by tracing 
the fate of a fraction of its most bound particles. In this 
way we can identify the point at which a halo dissolves as 
a result of disruption within a larger (sub)halo and form a 
branch consisting of its main progenitor at all earlier times. 
We discard such a branch if at the point at which this sub- 
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Progenitor mass function of Aquarius halos at z = 2 



10° 
IO" 1 

IO" 3 



IO" 4 
IO" 5 



aa a aq 



Sir 



o 



10" 6 
IO" 6 



This work O 

Aquarius CDM A 

WDM realization ■ 

WDM realization 1 □ 

WDM realization 2 □ 

WDM realization 3 ■ 



1 fl- 



lC" 4 10" 



10" 2 



io- 



M. 



progenitor/ ^^rinal 



10" 



Figure 6. Progenitor mass functions of four realizations of the Milky Way-mass WDM dark matter halos at z = 1 and 2 (squares) 
compared with the predictions of this work (circles). For reference, the corresponding progenitor mass functions from one realization of 
the Aquarius CDM simulations are shown as triangles. Each panel shows the fraction of the halo's mass contributed by progenitors in 
each mass bin. Masses are shown as a fraction of the final halo mass. 



halo had half its maximum mass the subhalo fails the cut 
on initial axis ratio. 

Fig. 6 compares the progenitor mass functions measured 
from a large number of merger tree realizations generated 
using the techniques described in this work (circles) com- 
pared with progenitor mass functions measured from four 
Milky Way-mass WDM simulation halos (squares). For ref- 
erence, the equivalent CDM progenitor mass function from a 
single realization of the Aquarius simulations is shown (tri- 
angles). Given the halo-to-halo scatter, there is good agree- 
ment down to Mp r0 genitor /Affi na i ~ 10 _ . Below this, the 
abundance of progenitors in the WDM N-body simulations 
exceeds that predicted by our techniques. At these mass 
scales, the artificial halo rejection algorithm reduces the 
number of progenitors by over an order of magnitude. This 
excess of low mass progenitor halos is consistent with fail- 
ure rate for the artificial halo rejection algorithm of slightly 
less than 10%. We cannot rule out such a failure rate in the 
algorithm and so the true progenitor mass function could 
decline much more rapidly. The progenitor mass functions 
from our analytic calculations decline rapidly as they ap- 
proach Mp r0 g en itor /Affinal ~ 10~ 4 , but then transition to a 
slower, power-law decline at lower masses. As will be dis- 
cussed in §3.3.3, this appears to be due to the finite differ- 
ence approximation used to compute merger rates (see §2.4) 
and so can be suppressed as necessary by lowering the value 
of e. 

In the CDM case the mass function (expressed 
in this way) levels off to be almost constant below 

Afprogenitor 

/M flnal ~ IO" 3 . In our calculation, which tracks 
the CDM case extremely well at high masses, this flatten- 
ing begins (for M pr0 gcnitor/Mfi na i > 10~ 3 ) but is then in- 
terrupted by the cut-off due to WDM physics. The result 
is a "bump" feature. The presence of such a bump is less 
clear in the WDM N-body progenitor mass functions (al- 
though there is a hint of it in the z = 1 results), due to the 
noisiness of those results. This feature may have interest- 
ing implications for the expected number of surviving dwarf 



scale subhalos, and so represents an interesting avenue for 
future investigation. 

3.2 Limitations of Only Modifying Power 
Spectrum 

Having demonstrated that our model is consistent with the 
available N-body simulation we now consider the effects of 
treating WDM incorrectly or incompletely in the extended 
Press-Schechter approach, as has previously happened in the 
literature (Menci et al. 2012). Fig. 7 shows a series of results 
for the halo mass function for 1.5 keV WDM. For reference, 
line 1 shows the mass function for CDM derived assum- 
ing a mass-independent critical overdensity. Line 2 switches 
to using the WDM power spectrum, but retains the CDM 
mass-independent critical overdensity. There is a clear sup- 
pression of the mass function below M B . Line 3 now uses the 
WDM power spectrum and adds mass dependence to the 
critical overdensity, but still uses the barrier crossing solu- 
tion for a linear barrier 11 given in eqn. (Bl). The suppression 
of the mass function is almost unchanged compared to the 
previous case where we used the CDM critical overdensity. 
Finally, line 4 shows the result of switching to using the cor- 
rect, numerically determined, functional form for f(S). The 
mass scale of the cut-off undergoes a dramatic shift of al- 
most an order of magnitude to higher mass compared to the 
previous case. The reason for this is simple to understand. 
When using the WDM B(S) in /linear (SO we are implicitly 
assuming that the barrier was flat at the same value of B(S) 
at all smaller 5* (this being the assumption used to derive 
/linear (5'). This gives a certain crossing probability. When 
we numerically determine /correct (S) the full S-dependence 
of the barrier is taken into account. In the first case, a ran- 
dom walk crossing B(S) at S must have remained below 

11 Specifically, we adopt a B(S) corresponding to the mass- 
dependent WDM critical overdensity, but simply use it in the 
solution for f(S) appropriate to a linear barrier. 
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Figure 7. The dark matter halo mass function for CDM (line 
1) and for various approximations to the 1.5 keV WDM solution. 
Line 2 shows the result of using the WDM transfer function, but 
retaining the CDM mass-independent critical overdensity and the 
functional form of f(S) appropriate to a linear barrier, /linear C^O- 
Line 3 improves on this approximation by using the correct mass- 
dependent WDM critical overdensity, but still uses fu nei3LY (S). Fi- 
nally, line 4 shows the result of using the correct functional form 
of f(S) for the WDM barrier. Note that, for large halo masses, 
all lines coincide precisely and so are hidden beneath Line 4. 



B(S) for all 5*' < S, while in the second case it must have 
remained below B(S') for all S' < S. Since B(S') < B{S) 
for all S' < S this second condition is much more strin- 
gent, and so far fewer random walks will satisfy it. Thus, 
the crossing probability, and so the mass function, will be 
more strongly suppressed. This illustrates the importance 
of correctly solving the barrier crossing problem for WDM 
calculations. 



3.3 Warm vs. Cold Dark Matter 

In the following we compare example results for CDM and 
1.5 keV WDM cases. This WDM particle mass differs from 
those used in previous sections (where the choice was con- 
strained to match those assumed in different N-body simu- 
lations). In all cases, we use the correct f(S) (determined 
numerically for the WDM case) and include remapping of 
the barrier for calculations of first crossing probabilities and 
mass functions. 



3.3.1 Mass Functions 

Fig. 8 shows the first crossing probability, f(S), for both 
CDM and WDM as a function of mass, M. The two are off- 
set by a constant multiplicative factor at large masses (small 
S). This is expected-we have chosen to re-scale the WDM 
barrier such that the halo mass function remains unchanged 
relative to the CDM expectation for large masses when we 
use the sharp fc-space filter in our WDM calculations. The 
mapping from f(S) to the mass function, n(M), is propor- 
tional to dS/dM: 




(20) 
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Figure 8. The first crossing probability, f(S), is shown as a func- 
tion of mass, M, for both CDM and WDM cases for 2 = 0. 



For large masses, dS/dM is larger at fixed mass for large 
masses in our WDM calculations compared to the equiva- 
lent CDM calculation (a consequence of the different win- 
dow functions adopted for the two cases). This difference 
in dS/dM offsets the difference in f(S) for large masses re- 
sulting in halo mass functions that agree between WDM and 
CDM. 

Below around M s the WDM first crossing probability is 
suppressed due to the rapidly rising barrier B(S). (There is 
a small region where the WDM f(S) exceeds that of CDM 
due to differences in the mapping from M to S in the two 
cases. 

These first crossing distributions translate directly to 
halo mass functions, as shown by Lines 1 and 4 of Fig. 7. As 
expected, the suppression in f(S) in WDM translates to a 
strong suppression in the mass function below around M s . 
At larger masses, the two are indistinguishable. 



3.3.2 Merger Rates 

Fig. 9 shows the rate of first crossing for CDM and WDM 
barriers for a 10 12 Mq halo at z = 0. The two lines almost 
coincide at high mass, although the WDM line lies slightly 
below the CDM line 12 . The merger rate in the WDM case 
is strongly suppressed below around M s due to the lack of 
halo in that mass range. The slight enhancement in the rate 
of first crossing in WDM compared to CDM just above M s 
is due to the different mapping between mass and variance 
in the two cases. 



3.3.3 Merger Tree Statistics 

Using merger rates computed as described in §2.4 we con- 
struct merger trees in both CDM and WDM cases using the 
algorithm described in §2.5, beginning with a halo of mass 
10 12 M© at z — 0. We generate 1743 trees in each case and 



12 For the same reason as in Fig. 8. Here the difference is less 
evident, partly because the two lines are closer to vertical, but 
also because the largest mass scale showed here is lower that in 
Fig. 8 such that S(M) differs less from WDM to CDM cases. 
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Figure 9. The rate of first crossing, df/dt(S), is shown as a 
function of mass, M, for both CDM and WDM cases, for the 
conditional barrier appropriate to a 10 12 Mq halo at z = 0. 

construct the mean progenitor mass function and the mean 
mass accretion history. 

Fig. 10 shows progenitor mass functions at z = 1, 3 
and 7. The sharp cut-off at 10 7 Mq (present also in the 
CDM case) is due to the imposed resolution of our merger 
trees 13 (and so is unphysical). With smooth accretion in- 
cluded, the WDM progenitor mass function closely matches 
that of CDM above about 3M a , but is strongly suppressed 
below it at lower masses. 

Well below the suppression scale a population of pro- 
genitor masses much less than M s builds up. These are the 
result of smooth accretion-the first crossing rate distribu- 
tion (see Fig. 8) cuts off below M B so there is no way for 
these halos to arise through branching of the merger tree- 
which gradually reduces the mass of the lowest mass halos 
going back in time. The numerical robustness of our model 
in this regime is discussed in Appendix B4, in which we also 
demonstrate that the position of the peak in the progenitor 
mass function is numerically robust and well determined. 

Fig. 11 shows the corresponding mean mass accretion 
histories of this halo (i.e. the mean mass of the most massive 
progenitor at each redshift). There is almost no difference 
between CDM and WDM. This is in agreement with the re- 
sults of Knebe et al. (2002) who found almost no difference 
in the mass accretion histories of individual halos in CDM 
and WDM N-body simulations. The two halos presented by 
Knebe et al. (2002) were significantly more massive (almost 
10 14 Mq) than those considered here, but they report the 
same conclusion for lower mass halos. Knebe et al. (2002) 
conclude that the number of mergers that contribute signifi- 
cant mass to the assembly of the halos is unchanged between 
CDM and WDM case. Our results suggest that this is an ac- 
curate conclusion for more massive halos (the assembly of 
which will be dominated by halos well above M s ). However, 
for lower mass halos, which gain a significant fraction of their 



13 The merger tree resolution is limited only by the available com- 
putational time and memory. We choose a resolution of 10 7 Mq 
in this case to be sufficiently below M s while keeping computing 
times tractable. 



mass from halos close to M a , our results clearly show that 
smooth accretion plays a crucial role in shaping the mass 
accretion history of WDM halos-without it, substantial dif- 
ferences from the CDM case would occur. 



4 DISCUSSION 

We have described algorithms for constructing halo mass 
functions and merger trees for dark matter halos in WDM 
universes. Our methods improve upon previous treatments 
which did not correctly solve the barrier first crossing prob- 
lem (Menci et al. 2012) and which used a top-hat filter to 
compute cr(M) resulting in an overestimate of the abun- 
dance of low-mass halos. Our results are in excellent agree- 
ment with the available N-body simulations. Illustrative re- 
sults clearly demonstrate that the mass function, and pro- 
genitor mass functions of WDM halos are strongly sup- 
pressed relative to CDM below about M s . Differences be- 
tween CDM and WDM in coarse-grained statistics, such as 
the mass accretion history, are small for halos well above 
the cut-off scale, providing that the accretion of smoothly 
distributed matter in WDM is accounted for. 

The method that we describe has a single free 
parameter — the coefficienct a appearing in eqn. (6). We have 
fixed the value of this parameter to match the location of the 
turnover in the N-body mass function reported by Schnei- 
der et al. (2012). This introduces a dependency on one of the 
simulations to which we compare our model. Nevertheless, 
as discussed in §2.1.1, we expect that the same value of a 
will be appropriate for all WDM particle masses of interest. 
The limited evidence available from our present work (i.e. 
that a chosen to fit the mass function for 0.25 keV WDM 
also successfully matches the progenitor mass functions for 
2.20 keV WDM) certainly supports this claim. 

Our approach has the advantage, compared to N-body 
simulations of WDM, of not being affected by numerical 
noise in the particle distribution which leads to the forma- 
tion of large numbers of artificial low-mass halos (Wang & 
White 2007), and of being substantially faster to evaluate 
mass functions and progenitor distributions. Using the tech- 
niques developed in this work, the procedure for applying 
them to dark matter with different phenomenology, or to 
other physics that modifies the power spectrum or excur- 
sion set barrier is straightforward: 

(i) Determine the linear theory power spectrum of density 
perturbations from the dark matter (or other) physics; 

(ii) Determine, through analytic calculation or idealized 
N-body simulations, the critical linear theory overdensity for 
collapse, which will depend on the (thermal and interaction) 
physics of the dark matter particle. 

Given these two inputs our techniques can be used to deter- 
mine the resulting halo mass function and merger histories 
of dark matter halos consistent with the input physics. The 
accuracy of our methods for phenomenology beyond that 
exhibited by WDM remains to be tested, but the success in 
this case leads us to expect that our methods will be gener- 
ally applicable. 

The nature of the dark matter distribution on small 
mass scales will be investigated by future lensing programs 
(Keeton & Moustakas 2009; Vegetti & Koopmans 2009a). 
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The techniques described in this work will allow detailed 
statistical predictions to be made for the expected numbers 
and masses of dark matter substructure as a function of dark 
matter particle properties. 

To make accurate predictions for the dark matter sub- 
halo distribution we have addressed only the first part of the 
problem, namely halo formation and merging. The second 
part, halo destruction by tidal forces must also be addressed. 
We plan to explore this process using the methods of Ben- 
son et al. (2002b; see also Taylor & Babul 2004), together 
with prescriptions for the internal structure of WDM ha- 
los (which will, of course, differ from that of CDM halos; 
Maccio' et al. 2012; Maccio et al. 2012). 

The methods described in §2 have been implemented 
within the open source semi-analytic galaxy formation 
code, GALACTICUS (Benson 2012). All results presented in 
this work were generated using GALACTICUS vO.9.1 r903. 
A control files and scripts to generate all results pre- 
sented in this paper using GALACTICUS can be found at 
http : / /users . obs . carnegiescience . edu/ abenson/ 
galacticus/parameters/dmMergingBeyondCDMl .tar .bz2. 
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Figure 11. Mean mass accretion histories derived via merger 
tree construction in CDM and WDM cases for a 10 12 M halo 
at z = 0. 



Figure 10. Progenitor mass functions derived via merger tree 
construction in CDM and WDM cases for a fO 12 Mq halo at 

2 = 0. 
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APPENDIX A: NUMERICAL METHOD 



In this section, we describe our numerical method for solving the excursion set barrier first crossing problem (eqn. 12). In the 
absence of a barrier, P(S, 5) would be equal to Po(S, 5) which is simply a Gaussian distribution with variance 5: 



Po(8,S) = 



1 



: exp 



25 



/2nS 

Since the barrier absorbs any random walks which cross is at smaller 5, the actual P(S, 5) must therefore be given by: 

rS 



P(6,S) = P (5,S)- f f(S')P [8-B(S'),S-S']dS'. 
Jo 



(Al) 



(A2) 



The second term on the right hand side of eqn. (A2) represents the distribution of random trajectories originating from the 
point (S,B(S)). The integral therefore gives the fraction of trajectories which crossed the barrier at 5 < S' and which can 
now be found at {S,8). 

Using this result, we can rewrite eqn. (12): 



1= I" f(S')dS'+ f B{S \p {8,S)- f f{S')P {5-B(S'),S-S')dS') 

JO J— oo L JO 



dS, 



in general and, for the Gaussian distribution of eqn. (Al): 



1 = / f(S')dS' + / 

JO J — oo 



/2tT5 



exp 



6^ 
25 



/ f(S') 
Jo 



The integral over dS can be carried out analytically to give 

" s r B(S) 



= / f(S')dS' 
Jo 



+ erf 



/(S')erf 



[5-B(S')Y 
y/2n(S-S') P V 2(5-5') 

B(5) - B(S') 



dS' 



dS. 



y/2(S-S') 



dS". 



We now discretize eqn. (A5). Specifically, we divide the 5 space into N intervals defined by the points: 

, ifi = 

S ' - 1 Er 1 A5 l if-' 1. 



(A3) 



(A4) 



(A5) 



(A6) 



Note that /(0) = by definition, so /(5o) = always. We choose ASi = S m nx/N (i.e. uniform spacing in 5) when computing 
first crossing distributions, and ASi oc 5; (i.e. uniform spacing in log(5)) when computing first crossing rates. 
Discrctizing the integrals in eqn. (A5) gives: 



f Sj f(S')dS' = £ 

J( > r=0 



f(Si) + f{S i+1 ) 



AS, 



and: 



F' /(S>rf 
Jo 



B(S) - B(S') 
y/2(S-S') 



2-1 



B(S 3 ) - B(Sj) 
^2(5, - 50 



+ /(5 i+1 )erf 



5(5,0 ~ B{Si + i) 
y/2(Sj - 5 I+1 ) 



ASi 



We can now rewrite eqn. (A5) in discretized form: 



3-1 



f(Si) + f(S i+1 ) 



ASi + erf 



B(Sj) 
\/257 



3-1 



B{Sj)-B(Sj 
y/2(Sj - Si) 



+ /(5 i+1 )erf 



B(Sj) - B(S i+1 ) 
^2(5, - Si+i) 



Solving eqn. (A9) for f{Sj): 



(A7) 
(A8) 

ASi. 
(A9) 



-erf 



B(Sj) - B(Si) 



3-2 



A5 J -i/(5 J ) 



1 _ j2 f(Si)+m+i) ASi _ nsj-i) AS ._ i _ erf 



i=0 
3-2 



2 2 
[BiS,) - B(Si)] 



+ 



2 3 \ \/2(^-5,_0 / 



+ /(5 i+ i)erf 



B(Sj) 
\/2S~i 



[B(Sj) - B(S i+1 ) 
\/2(5j — Si+i) 



A5 



3-1- 



ASi 



(A10) 



For all barriers that we consider: 



erf 



B(S 3 ) -B(S 3 ) 



(All) 
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We can then simplify eqn. (A10) 

2 f(s.) ■ - • 

i=0 



f(Si) = 



AS,- 



'3-1/ 



AS,_i - erf 



,)-B(S i+ i)] 
\/2(S 3 - Si+i) 



+ 2/(^_i)erf 



[5(5,-) - 5(5,-1)] 



\/2(5j - Sj-i) 
Consolidating terms in the summations: 



AS 



3-1 



AS 



3-1 



1 - erf 



ga-erf 



In the case of constant ASi(= AS) this can be simplified further: 



f{Sj) 



AS 



1 - erf 



B(Sj) 
\/2S7 



2 ga-erf 



ASi 



f(Si 



BjSj) - B(Sj) 
V 2 (Sj ~ Si) 



AS^ + AS, 



f(Si). 



(A12) 



(A13) 



(A14) 



In either case (i.e. eqns. A13 and A14) solution proceeds recursively: /(So) = by definition, /(Si) depends only on the 
known barrier and /(So), f(Sj) depends only on the known barrier and /(S<j). 
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Figure Bl. The absolute fractional error in the excursion set 
first crossing distribution function, f(S), computed using our nu- 
merical method compared to the analytic solution for a constant 
barrier. 



APPENDIX B: NUMERICAL TESTS 

In this Appendix we examine the numerical accuracy and 
robustness of our methods. 



Bl First Crossing Probability Solutions 

To test the accuracy of our numerical solver for the first 
crossing distribution we compare its results to the known 
analytic solution for a constant barrier. Specifically, we con- 
sider the CDM case with no remapping of the barrier such 
that B(S) — S c ^constant and 



f(S) = 



exp 



5 ± 

2S 



(Bl) 



SV2tyS 

Fig. Bl shows the fractional difference between results ob- 
tained using our numerical method and the analytic solution 
in this case. The numerical error is small except for at the 
highest masses (smallest variances) where the discreteness 
in our grid becomes an issue and the error reaches a few 
percent. The accuracy achieved is sufficient for the present 
work (where we are mostly concerned with low-mass sys- 
tems) and is easily improved by adopting a smaller value of 
AS. 

B2 First Crossing Rate Solutions 

Once again, to test the accuracy of our numerical solver 
in the case of computing first crossing rates we compare its 
results to the analytic solution for a constant barrier. Adopt- 
ing the same constant barrier model as in Appendix Bl the 
crossing rate is 

df _ 1 1 dS c 

dt ~ (S - So) 3 / 2 d*~ 

where So is the variance corresponding to the mass of the 
final halo (i.e. the halo formed through the merger event). 
Fig. B2 shows the fractional error in the numerically derived 
merger rate compared to this analytic solution for e = 0.01. 
The fractional error is constant across the range of masses 
shown and is limited by the value of e chosen. 



(B2) 
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Figure B2. The absolute fractional error in the excursion set first 
crossing rate function, df(S)/dt, computed using our numerical 
method compared to the analytic solution for a constant barrier. 



B3 Testing the Parkinson-Cole-Helly Algorithm 
on Smaller Mass Scales 

The Parkinson et al. (2008) empirical modification to the 
merger tree branching rate was calibrated against N-body 
merger trees drawn from the Millennium Simulation. As 
such, it has been tested for masses above approximately 
1O 1O M0. Here we employ this same modification for much 
lower masses. Fig. B3 compares progenitor mass functions 
generated by the Parkinson et al. (2008) empirical modifica- 
tion with those extracted from the Aquarius CDM simula- 
tions of Springel et al. (2008) which resolve halos of masses 
1O 6 M . It can be seen that the Parkinson et al. (2008) em- 
pirical modification performs very well (with unchanged pa- 
rameter values) for these much lower mass halos in the CDM 
case also. 



B4 Merger Tree Accuracy and Smooth Accretion 

To test the convergence of our merger trees with respect 
to the parameter e used in our numerical determination of 
the first crossing rate distribution (see §2.4) we compute 
progenitor mass functions of 10 12 Mq halos for e = 0.010, 
0.003, and 0.001. Additionally, we perform these calculations 
both with and without the smooth accretion term of §2.5.1. 

Fig. B4 shows progenitor mass functions at z = 1, 3 and 
7. When smooth accretion is included, the WDM progenitor 
mass function closely matches that of CDM above about 
3M S , but is strongly suppressed below it at lower masses. 
If smooth accretion is ignored the WDM progenitor mass 
function evolves much more slowly and diverges from the 
CDM case even at the highest masses. Ignoring this smooth 
accretion leads to significantly biased results. 

It can clearly be seen that our WDM results are con- 
verged with respect to the e parameter except at masses well 
below the suppression scale in the progenitor mass function 
in cases where smooth accretion is included. Here, a pop- 
ulation of progenitor masses much less than M s builds up. 
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Progenitor mass function of CDM Aquarius halos at z = 1 
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Progenitor mass function of CDM Aquarius halos at z = 2 
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Figure B3. Progenitor mass functions of the Aquarius CDM dark matter halos (Springel et al. 2008) at z = 1 and 2 (histograms) 
compared with the predictions of the Parkinson et al. (2008) algorithm. Each panel shows the fraction of the halo's mass contributed by 
progenitors in each mass bin. Masses are shown as a fraction of the final halo mass. Error bars on the predictions from the Parkinson 
et al. (2008) algorithm indicate the 20 th and 80 th percentiles of the distribution of progenitor mass functions. While this algorithm was 
calibrated on the much higher mass halos found in the Millennium Simulation (Springel et al. 2005) it can be seen to work extremely 
well at these lower masses also. 



These are the result of smooth accretion-the first crossing 
rate distribution (see Fig. 8) cuts off below M s so there is 
no way for these halos to arise through branching of the 
merger tree-which gradually reduces the mass of the lowest 
mass halos going back in time. 

Our numerical determination of the first crossing rate 
function is currently not robust in its determination of 
smooth accretion rates for these lowest mass (highest vari- 
ance) halos where the excursion set barrier is a very rapidly 
changing function of variance and our discretization of S 
used to obtain a numerical solution inevitably does a poor 
job of resolving the barrier. This could of course be mitigated 
by using a yet smaller value of e (requiring substantially in- 
creased precision in the numerical solutions) or a finer grid in 
5* (requiring both increased precision and substantially more 
computing time). Nevertheless, the position of the peak in 
the progenitor mass function is well determined, and the 
differences in the abundance of low mass progenitors have 
negligible effect on the mean mass accretion histories of ha- 
los considered in this work. 
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Dark Matter Halo Merger Histories Beyond Cold Dark Matter 
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Figure B4. Progenitor mass functions derived via merger tree 
construction in CDM and WDM cases for a fO 12 Mq halo at 
z = 0. For the WDM case, results are shown for different values 
of e and also for cases with and without the accretion of smoothly 
distributed matter. 
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